##################################################################################
##################################################################################
####Non-state Atrocities in Capital Cities - War Type and Insurgent Atrocities####
##################################################################################
##################################################################################
library(MASS) 
library(pscl)
library(foreign)
library(stargazer)
library(mvtnorm)
library(plyr)
library(interplot)


# Set working library
setwd("~/Data/Conflict Type Analysis/")
# Read in main data
war.etsec.dat <- read.dta("eth_sec_war_dat.dta")

#########################################
###ZINB Model Corresponding to Model 4###
#########################################

####Baseline model without interaction
at.zinb.4.a <- zeroinfl(incidentnonstatefull ~ lag_Capital +
                          lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                        + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                          year98 + year99 + year00 + year01 + year02 + year03 + year04
                        |lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                        + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea, 
                        dist = "negbin", data = war.etsec.dat)
summary(at.zinb.4.a)
AIC(at.zinb.4.a)

####Ethnic wars
at.zinb.4.ce <- zeroinfl(incidentnonstatefull ~ lag_Capital + ETHNOWAR + lag_Capital*ETHNOWAR + 
                           lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                         + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                           year98 + year99 + year00 + year01 + year02 + year03 + year04
                         |lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                         + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea, 
                         dist = "negbin", data = war.etsec.dat)
summary(at.zinb.4.ce)
AIC(at.zinb.4.ce)

####Wars of secession
at.zinb.4.cs <- zeroinfl(incidentnonstatefull ~ lag_Capital + SECESSION + lag_Capital*SECESSION +
                           lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                         + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                           year98 + year99 + year00 + year01 + year02 + year03 + year04
                         |lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                         + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea, 
                         dist = "negbin", data = war.etsec.dat)
summary(at.zinb.4.cs)
AIC(at.zinb.4.cs)


##############################################################################################
###Logit Models For Predicted Impact Of War Type On Insurgent Propensity To Target Capitals###
##############################################################################################
####Create a binary variable atrocity/no atrocity
war.etsec.dat$mass_inc <- ifelse(war.etsec.dat$incidentnonstatefull>0,1,0)

####Baseline model without interaction
at.log.4.a <- glm(mass_inc ~ lag_Capital +
                     lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                   + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                     year98 + year99 + year00 + year01 + year02 + year03 + year04,
                   family="binomial", data = war.etsec.dat)
summary(at.log.4.a)
AIC(at.log.4.a)

####Ethnic wars
at.log.4.ce <- glm(mass_inc ~ lag_Capital + ETHNOWAR + lag_Capital*ETHNOWAR +
                     lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                   + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                     year98 + year99 + year00 + year01 + year02 + year03 + year04,
                   family="binomial", data = war.etsec.dat)
summary(at.log.4.ce)
AIC(at.log.4.ce)
##Plot the conditional effect of ethnic war on capital city's coefficient
pdf("intce.pdf")
interplot(m = at.log.4.ce, var1 = "lag_Capital", var2 = "ETHNOWAR")+
  # Add labels for X and Y axes
  xlab("Ethnic War") +
  ylab("Estimated Coefficient for\nCapital") +
  # Change the background
  theme_bw() +
  # Add the title
  ggtitle("Estimated Effect of Ethnic War \non Capital's Impact on the Frequency of Insurgent Atrocities") +
  theme(plot.title = element_text(face="bold")) +
  # Add a horizontal line at y = 0
  geom_hline(yintercept = 0, linetype = "dashed")
dev.off()

####Wars of secession
at.log.4.cs <- glm(mass_inc ~ lag_Capital + SECESSION + lag_Capital*SECESSION +
                     lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                   + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                     year98 + year99 + year00 + year01 + year02 + year03 + year04,
                   family="binomial", data = war.etsec.dat)
summary(at.log.4.cs)
AIC(at.log.4.cs)
##Plot the conditional effect of secessionist war on capital city's coefficient
pdf("intcs.pdf")
interplot(m = at.log.4.cs, var1 = "lag_Capital", var2 = "SECESSION")+
  # Add labels for X and Y axes
  xlab("Secessionist War") +
  ylab("Estimated Coefficient for\nCapital") +
  # Change the background
  theme_bw() +
  # Add the title
  ggtitle("Estimated Effect of Secessionist War \non Capital's Impact on the Frequency of Insurgent Atrocities") +
  theme(plot.title = element_text(face="bold")) +
  # Add a horizontal line at y = 0
  geom_hline(yintercept = 0, linetype = "dashed")
dev.off()

####Export all models to LaTex
stargazer(at.zinb.4.a, at.zinb.4.ce, at.zinb.4.cs, at.log.4.a, at.log.4.ce, at.log.4.cs)
stargazer(at.zinb.4.a, at.zinb.4.ce, at.zinb.4.cs, zero.component = T)
AIC(at.zinb.4.a, at.zinb.4.ce, at.zinb.4.cs, at.log.4.a, at.log.4.ce, at.log.4.cs)